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Abstract 



In shock precursors populated by accelerated cosmic rays (CR), the CR return 
current instability is believed to significantly enhance the pre-shock perturbations 
of magnetic field. We have obtained fully-nonlinear exact ideal MHD solutions 
supported by the CR return current. The solutions occur as localized spikes of cir- 
cularly polarized Alfven envelopes (solitons, or breathers). As the conventional 
(undriven) solitons, the obtained magnetic spikes propagate at a speed C propor- 
tional to their amplitude, C = C^max/ V^Bq. The sufficiently strong solitons run 
thus ahead of the main shock and stand in the precursor, being supported by the re- 
turn current. This property of the nonlinear solutions is strikingly different from 
the linear theory that predicts non-propagating (that is, convected downstream) 
circularly polarized waves. The nonlinear solutions may come either in isolated 
pulses (solitons) or in soliton-trains (cnoidal waves). The morphological similar- 
ity of such quasi-periodic soliton chains with recently observed X-ray stripes in 
Tycho supernova remnant (SNR) is briefly discussed. The magnetic field amplifi- 
cation determined by the suggested saturation process is obtained as a function of 
decreasing SNR blast wave velocity during its evolution from the ejecta-dominated 
to the Sedov-Taylor stage. 
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Subject headings: acceleration of particles — cosmic rays — shock waves — ISM: 
supernova remnants — magnetohydrodynamics (MHD) — magnetic fields 



Introduction 



The nonresonant cosmic ray (CR) return current instability (also termed as Bell's instability) 
is expected to bootstrap the acceleration of CR in shocks by enhancing the magnetic field in 
the shock precursor. The most unstable is a circularly polarized, field alig ned, aperiodic mode , 



similar to the internal kink (Kruskal-Shafranov) mode in plasmas (see, e.g., 



Rvu tov et al. 



In the contex t of th e CR acceleration in shocks, it was also studied by lAchterbergl (LL983); 



Shapiro et al 



(1998 



2006). 



Belli (12004)) reawakened the interest in this instability by emphasizing its role in magnetic 
field amplification and suggested its saturation due to magnetic tension. But, since the growth 
rate decreases with the wave number only as Vk, the magnetic tension term does not stabilize 
the long waves. This opens the door for a strong, 8B ^> Bq field amplification. The caveat is that 
the non-propagating long waves have limited (precursor-crossing) time to grow. By contrast, 
the nonlinear solutions, that we present in this paper, can stand off in the flow ahead of the 
shock, thus warranting the saturation. 



Vladimirov et al. 



There have been considerable efforts to understand the Bel 
with a strong emphasis on t he MH D and PIC simulations 



2006; 



Reville et al. 



2007; 



Niemiec et al 



2008 



's mode saturation mechanisms, 



(Bell 



2004 



Pelletier et al 1 12 006 



Zirakashvili et al. 



2008 



Bvk ov et al. 



2009 



Luo & Melro se 2009; 



Riquelme & Spitkovskvl2009 : 



Stroman et al. 



Bell (2005) (see also 



2009 



Dieckmann et al. 



Niemiec et al 



2010). As first demonstrated in 3D MHD simulations by 
2008), the saturation is achieved when the Ampere force expels plasma and the helical mag- 
netic field radially, thus forming plasma cavities. The instability appears to saturate only in 
3D, or at least requires a quasi-2D dynamics, perpendicular to the ambient field. However, the 
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fastest growing modes are field-aligned, i.e., at least initially one-dimensional. Therefore, it 
is necessary to understand structures that form at the ID phase and particularly the nonlinear 
mechanisms of their saturation and propagation ahead of the shock. These structures may in the 
main cease to grow before the subsequent 3D dynamics kick in largely by spreading the satu- 
rated turbulence energy in £-space. Although this scenario may appear to be at odds with many 
simulations, recent Chandra observations of the Tycho supernova remnant (SNR), for exam- 
ple, indicate the presence of quasi- ID structures (s tripes), inconsisten t with the quasi-isotropic 



201 Moreover, while being 



nonlinear dynamics observed in those simulations (|Eriksen et al. 
very useful for our understanding of CR instabilities, simulations cover only a tiny fraction of 
the dynamical range of typical SNR-shock acceleration process and introduce artificial dissipa- 
tion in essentially collisionless plasmas. 

Alfven waves usually saturate by modulational instability. However, being a strong MHD 
aperiodic instability, Bell's ins tability hampers direct ap plications of standard methods, such 



as the weak-turbulence theory (|Sagdeev & Galeev 



19691). The latter typically deals with prop- 



agating and weakly interacting eigen modes and, as a driver amplifies them, they cascade the 
wave energy to the dissipation scale. The Bell's linear mode does not propagate (in the linear 
approximation), and does not even exists without the driving current. The lack of long wave 
stabilization is also based on the comparison of linear contributions to the square of the growth 
rate of the driving current (oc k) and magnetic tension (°c k 2 , eq. |[T2l below). A clue to saturation 
in a similar syste m of the pressure-a n isotro py-driven fire-hose instability is provided by an exact 



solution due to 



Berezin & Sagdeev 



(|1969|) . While at peaks of magnetic energy it takes nearly 
all the instability free energy (B\ /Sn ~ Pn — P±_ 3> Bq/Sji), on the average only the moderate 
field amplification ~ Bq is observed. 

In this paper we present an exact solution of the current-driven MHD equations (e.g., CR re- 
turn current). It differs from the linearly growing solution in that it propagates with the velocity 
proportional to its (constant) amplitude and is spatially localized. 
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2. Basic equations 



The linear theory of Bell's inst ability indi cates that the fastest growing modes are directed 



along the ambient magnetic field (IBell 



2005). Therefore, we will consider ID magnetohydro- 



dynamic (MHD) equations in a CR shock precursor using a coordinate system with the axis x 
along the field. The general ID equations read: 



dt ^ ^ ' dx x 



(1) 



dU x __ Id B\ 
dt pdxSn' 



(2) 



dV L B dB ± 1 

— ; — = z =; — I J x K - 

dt Anp ox cp 



(3) 



dB 



Br 



B 



dU x 



(4) 



dt " dx dx 
Here dj dt = d /dt + U x d/dx, p is the gas density, U x ,Bq and U^,B^ are the gas velocity and 
magnetic field components along the field and in the (y, z)-plane, respectively. The x-component 
of magnetic field B x = Bq = const because of V ■ B = 0. In eq.© we have included the plasma 
return current by representing the total plasma current as J tot = (c/Atz) V x B + J, where the 
part of the plasma current J = —Jcr compensates the CR current. Eq.© implies that in our 



reference frame J x Bq = 0. We neglect the thermal and CR pressure, as 



Bell 



(2004) did. It 



should be noted, however, that these CR pressur e gradient drives an aco ustic instability of the 



shock precursor (also called Drury's instability, 



Drury 



19841 : 



Dorfi 



1984). Moreo ver, the acous 



Malkovetal 



tic i nstability grows faster than the Bell's instability for /3 = SkP/Bq < 1 (see 
2010 where the studies of evolution, saturation, as well as the associated particle transport and 
cascading of magnetic energy, are also referenced). 
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Returning to eqs.([TH4]), it is convenient to introduce a Lagrangian mass coordinate t, : 

d^ = ^- (dx - U x dt) , (5) 
Po 

where po is the background density. Considering scales shorter than the precursor size, we treat 
Po and the bulk plasma speed UxQ as coordinate independent (U x q = in the plasma frame). 

Next, we reduce eqs.£[]|4]) to the following system of two equations that describe the magnetic 
field and density perturbations: 



d 2 B 2 d 2 B i d B 

df 2 p- CA W7o ~ ^o BoJ Mp (6) 

d 2 Po 2 d 2 \B\ 2 

— H — — = (7) 

where 



B = B y + iB z and C% 



2 _ B 



4^p 

The r.h.s. of eqs.© is the instability driver. Without it, the equations describe the conventional 
MHD modes, propagating at an arbitrary angle to the ambient magnetic field. By choosing the 
averaged components B y = B z = 0, we restrict our treatment to the parallel propagation along 
the x-direction. 

3. Traveling wave solutions 

We look for the solutions of the system given by eqs.© and CO) in the form of a traveling 
wave: 

B = B max v(Qe- icot (8) 
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p = p(C) 

where £ = t, — Ct, C is the (constant) propagation speed of the traveling wave, 5 max is the wave 
amplitude that we specify in eq.© below, and co = 9t ft) is the wave frequency. Note that for 
ft) 7^ 0, the solution is not steady in any reference frame. For that reason, the spatially localized 
version of this solution is some times called 'breather' as opposed to the soliton, customary to 
CO = case. Integrating then eq.© twice, we obtain 



"° l-l£ (9) 



max 



P Bl 

where 5 max = SnpoC 2 . We have chosen the integration constants in such a way that B — > for 
p — > po (background plasma) and B — > B max for p — > °° (flow stagnation point, if present). This 
sets the interval for variation of v (£) : < v < 1. 

Substituting B from eq.® and p from the last equation, eq.© yields 



v| 2 ) v - iK^r fl " M 2 ) v - ^ f 1 - \v\ 2 ) V = 0. (10, 



d£ 2 V 1 1 ) dC V 1 1 / C 2 

Here we have used the following notations 



BrJ CO Bl C\ 

^=^7^- 2 7T, a= 1-2-^- = 1--4, (11) 
cp C 2 C 5 2 fl * C 2 ' 

where = B^/AitpQ. The linear dispersion relation can be recovered by letting v(£) °= 
v -> in eq.([T0l): 



co = kC±^k 2 Cl+B Jk/cp . (12) 

The arbitrary propagation speed C is a parameter of a Galilean transformation (zero in the 
plasma frame), while the imaginary part of ft) is an invariant of such transformation as it should 
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be. In the nonlinear treatment the wave velocity with respect to the plasma depends on the wave 
amplitude (nonlinear dispersion relation). Meanwhile, the linear instability occurs in the long 
wave limit for BqJIc < 0. It should be emphasized that only if the quadratic B term is neglected 
in eq.©, is there no coupling to the density modulations in eq.©. It is interesting to note that 
in the strong nonlinear limit B max /Bo — > °°, eq.dlOl) degenerates into a linear equation for the 
function v ^1 — |v| 2 j . This limit, however, cannot be understood without the nonlinear solution. 

To find such solution, we write 



v(Q = V^e i@ (13) 

where w(Q > 0. Substituting v from eq.([T3l into eq.dTOl) and separating the imaginary part, for 
the phase © we obtain the following equation: 

d® wP(w)+A 
as w(a — w) z 

where 

P = w 2 -(3a + l)w/2 + a 

We have introduced a new variable s = KC, /2 and an integration constant A. We may choose it 
by specifying the properties of the solution sought. The regularity of at w = implies A = 0. 
Next, taking the real part of eq.dTOl) and using eq.(fT4l) with A = 0, for w (s) we obtain 



2 n2 



d L f w z P 



ds 2 f 3 



+ 2^(1 



w) 



WP CD 1 







(15) 



where we have denoted 



f(w) = \/w(a — w) 
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Eq.(fT51) can be readily integrated by multiplying it by df/ds. We choose the integration constant 
to select an isolated pulse (soliton) solution of eq.([T5T), i.e. w — > 0, as s — > ±°°. Then, the first 
integral reads 



\as j (3w — a) (a — w) „=o 

where 

C = 4a 2 (ail 2 - l) 

C\ = 2a [2 (3+ a) -a/l 2 (7 + a)] 

C 2 = $H 2 a(a + 2)-a 2 -Ha-9 

C 3 = 2[2{3 + a)-jl 2 {3 + 5a)\ 
C 4 = 4(^2-1) 



with 



2^2 



[l l = A(o l /K l C 



2cp Cco ) 



Sagdeev 



1966) sug- 



A useful analogy between nonlinear waves and nonlinear oscillators (e.g., 
gests to interpret the first term eq.(fT6l) as kinetic and the second term as potential energy. The 
'oscillator's coordinate' w > 0, as a function of 'time' s, leaves w = 0+ at s = — °° and re- 
turns there at s = +°°. Periodic solutions (cnoidal waves) can also be easily by changing the 
integration constant. 

The amplitude wo(a,fi) of the localized solution (soliton) is obviously determined by the 
smallest positive root of the polynomial in eq.(fT6l). so that the "oscillator" bounces between 
w = and w = wq. In the simplest case of a small amplitude solution 



wo ~ — Cq/C\ a < 1 
(where Cq > and C\ < 0) the solution has a classical soliton profile 



Apart from the condition a\l 2 > 1 (to ensure wq > 0), i.e. 



1-CJ/C2< 



BqJ 



<l + v /l-C2/C 2 



IcpQCoC 



and the technical restriction a/i 2 — 1 <C 1 (to neglect the n > 1 terms in eq. |[T6l ). this solution 
imposes no further constraints on co and C. However, it has a very strong amplitude limitation, 
ajl 2 -1 < 1 (virtually a wave packet of linear waves). We are interested in an opposite case 
of highly superalfvenic solitons that are not convected rapidly with the flow into the sub- shock 
and may stay ahead of it, when C ^> Ca- 

A relation between CO and C (nonlinear dispersion relation) arises from the extension of the 
above solution to larger wq. Clearly, we have to pass the point w = a/3 smoothly which requires 
a double root of the polynomial in eq.(fT6l) at w = a/3: 



Interestingly, the both conditions are met simultaneously as soon as the following dispersion 
relation is satisfied: 



Recalling that the small amplitude soliton wq <^ 1 branches off from the trivial solution at the 
threshold a/i 2 = 1, in the case of Ca C, i.e. a ~ 1, we may accept eq.([T8l) to be valid in the 
entire parameter range afl 2 > 1 . 

Let us rewrite the above dispersion relation as follows 




k,C 



oj = 
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where we have defined the linear instability wave number (see eq. lfT2lO as kj = ItzJ /cBq. 
Strong solitons with Ma = C/Ca S> 1 have either high or low frequency: CO = (16/9) kjC, 
CO = kjC/2M\. The spatial scale of the solitons, given by the 'wave number' K, eq. (fTT|) . can be 
represented as follows, Fig{T] 



v { 2kj V(l-M A - 2 )/(l + l/8Mj) 
M A1±,/(1-M A 2 )/(1 + 1/8M2) 



It is interesting to note that both solutions disappear (spread to infinity) in the limit J — > 0, 
although they have disparate scales, particularly for Ma ^> 1. Therefore, the external current is 
essential and there is no transition to conventional simple wave MHD solutions for the vanishing 
CR-current. 

To obtain the spatial structure of the above solutions, we substitute eq.(fT8l into eq.([T6l). The 
latter takes the following simple form 



( — J =^r— z Z2Q (w) (20) 



\ds J 2a (a 



w) 



where 



Q 2 = w 2 -2hw + a; h = (a + 3)/4. 
Eq.(|20l) can be reduced to a quadrature: 



2a 

slw) 



I — a 

where 



, 1 h — w /-^ 
cosh — == + vai? 

yh 2 — a 



(21) 



_ , + , y/a — y/h 2 — a + h 
R = ln—= In . 

^a-w + Q y/a + y/h 2 -a-h 

Using eqs. (fl4)) and (T20l) . the solution for the phase (w) can be reduced to another quadrature 
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2a 
1 — a 



cosh 



i h — W „ . r- 



tan 



-l 



Q — w + a 



tan 



-l 



y/h 2 — a — h + a 



(22) 



'a(l -a) /I) \ y/a(l -a) /2 

The 5 x -component of the solitary solution is shown in Figj2](the e~ m1 -factor omitted). The 
wave packet in the compressed area becomes more oscillatory, as may also be seen from FigfJl 
which shows the soliton phase as a function of dimensionless coordinate s. The local dimen- 
sionless wave number stays approximately constant (d&/ds ~ 1), apart from the above phase 
steepening near the maximum amplitude, where d&/ds ~ 2. 



4. Maximum Magnetic Field 

The isolated solitons obtained in this paper belong to a one parameter family; the ampli- 
tude B max or Mach number Ma = B max / \/2Bq can be used as such parameter. In a CR shock 
precursor the soliton scale is determined by the scale of seed waves for the subsequent non- 
resonant instability. The seed waves are resonantly excited upstream of the strong CR current 
zone by the high energy CRs. Then, K ~ rj (/>*), where r g is the gyroradius of the seed 
generating CRs of momentum p = p*. This amounts to M\ = kjr g (p*) for the upper (long- 
wave) soliton branch in eq. (|T9b. Note that may be <C p max due to a poor CR confinement in 



the range p* < p < p max (|Malkov & Diamond 



2006). If the CR current is sufficiently strong, 



kj r g(p*) ^ 1 an d only the upper-sign soliton in eq.(fT9l) and Figfflcan accommodate the re- 
quirement K ~ r g 1 (p*) for Ma 3> 1 . Then, the maximum magnetic field for a given soliton can 
be written as B^^/Bq = 2M\ w 2kjr g (/?*), or 



5 max = *nV s ncRP*, 



(23) 
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where hqr is the CR density. 



The scaling of CR-enhance d magnetic 



energy with the ambient density p and shock velo city 



Belli (|2004 suggested B 2 /p « V s 3 , whereas 



Volk et al. 



(20051) in 



V s is debated in the literature, 
dicate that B 2 / p V 2 . The difference between the two scalings is whether a fraction of mechan- 
ical energy flux or momentum flux goes into magnetic e nergy . By co ntrast, eq.(l23l constitutes 



Vink (2008) summarizes the informa- 



the conversion of CR energy flux into magnetic energy, 
tion about the magnetic field from a number of SNR, with the two phenomenological scalings 
superimposed, Fig HI Note that 5 max in eq. ((2~3l) coincides with the condition of magnetization 
(trapping by the wave) of the current-carrying particles kj (B max ) r g (/>*,i? max ) = 1 /2, which is 
also (formally) similar to the Bell's phenomenological condition of balancing the Ampere force 
and the magnetic tension for the instability saturation. However, the saturation mechanism be- 



hind eq. (|23l) is different in that B max is only the peak mag netic field. The magnetic energy 



Berezin & Sagdeev 



1969). More 



density would be smaller by a soliton filling factor f s (cf. 
importantly, the efficiency of CR acceleration and subsequent conversion of their energy into 
magnetic field should depend on 7, V s and other acceleration parameters which almost certainly 
rules out the single power-law relation between B 2 /p and V s . 

Therefore, we obtain such relation in a different way, which we outline below and will de- 
scribe in detail elsewhere. Consider a nominal SNR with the shock speed V s slowing down from 
an initial V s =Vq= 1.34 x Vst* to V s ~ 0.1 x Vh (i.e., well into th e Sedov-Taylor phase) where 



V S T = 10400 x e\( 2 (M e /M Q ) l/z km/s dMcKee & Truelove 



1/2 , 



1995). Here E51 is the explosion 



energy in I0 5l ergs and M e -the ejecta mass. During its evolution, the SNR should follow the 
points sampled from a set of sup posedly similar remnants in FiglU Using the V s (t) dependence 



from (McKee & Truelove 



regime from (jMalkov & Drury 



1995), the momentum p* ~ p m ax(0 in the nonlinear acceleration 



200 lL eq.r7.45 



Next, we obtain hcr from eq.(15) in (|Malkov 



), we express p* in eq.(l23l as a function of V s 



1997|) for the evolving subshock strength with 



the particle injection rate held approximately constant in the efficient acceleration regime (see 
eq.(37) in the same reference). Using these results, we finally obtain from eq.(l23l an expression 
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for 5max/P' a g am > as a function of V s . A preliminary example of such calculation is shown in 
Fig J4] with the green line. In an intermediate range of V s the scaling is B^/ p V s 4 (close to 
the Bell's scaling) but it rolls over to turn to zero at V s = Vq ~ 9 • 10 4 km/ s. This is because the 
magnetic field generation is pinned to the CRs req. (|23l) l which are not yet there at V s = Vq, i.e. at 
t = 0. The other strong deviation from a power-law should occur at lower V s where acceleration 
bifurcates into a inefficient (test particle regime) through a characteristic S-curve. 



5. Discussion 

The purpose of this paper was to understand the nonlinear evolution of the non-resonant cur- 
rent driven instability by studying saturated nonlinear waves (solitons) as their ensemble (or 
that of their shock counterparts, if dissipation is efficient) may comprise the as ymptotic state of 



the sy stem. Such scenario is supported by simpler (but fully integrable, e.g., iKaup & Newell 



19781) weakly no nlinear MHD models, s uch as the derivative nonlinear Schroedinger equation 



(DNLS, see also 



Mjolhus & Hada 



1997 



for a review). In such models, arbitrary initial condi- 
tions evolve into an asymptotic state of quasi-independently propagating solitons, very much 
similar to those found in the present paper. The difference, however is that our system is driven 
by the CR return current and its solutions do not transition into the MHD solutions. 

The relevant question of soliton stability should be addressed in 2- 3D setting. The 2-3D 



2003) 



instability of a ID soliton could comprise a wave front self-focusing (jPassot & Sulem 
and thus elucidate the subsequent 3D structures. Such studies are beyond the scope of this letter, 
but a qualitative stability examination is in order. It may be based on the nonlinear dispersion 
relation given by eq.([T9l and FigfTJ The parts of the dispersion curves with d \K\ jdC < 
(where K and C °= Ma are the wave number and propagation speed) correspond to the solitons 
with negative dispersion and should be stable. The oft-used justification of stability is that a 
nonlinear steepening of the soliton's leading edge generates higher wave number modes and 
they should not run faster than the soliton itself (negative dispersion is thus required). It should 
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be also noted here that, once the soliton solutions of the driven system are obtained, they can 
also be arranged in a quasi-periodic or even chaotic soliton lattice. By add ing weak diss ipation, 



the leading edges of these solitons can be converted into shock fronts (ISagdeevlll966|) which 



usually increases the dissipation of the driver energy, thus reducing the saturation amplitude. 

To conclude, there exists upper bound on 5 max since solitons with C/Ca = B max / y/lBo > 
V s /Ca outrun the shock and cannot be sustained by the return current. However, as transients, 
they may promote particle acceleration far upstream to synergistically supply themselves with 
the CR return c urrent. This might be a plausible scenario for m uch-discussed CR acceleration 



bootstrap (e.g., 



Mal kov & Drurv 



2001 



Blandford & Funk 



2007b . Furthermore, strong solitons 



running ahead of the shock may becom e visible in X-rays as quasi-periodic stripes, similar 



to those recently observed by Chandra (|Eriksen et al. 



201 1|) in some parts of the Tycho SNR. 



The Eriksen's identification of the stripe spacing with the maximum gyroradius of acceler- 
ated particles is consistent with our determination of the soliton spacing in Secj4] but with a 
lower than 2 PeV energy. The scale is set by the highest energy particles ahead of the field 
amplification zone. The soliton wave length should be noticeably shorter than the distance be- 



tween them. At the same time, similar structures may resu 
the CR-pres sure-driven acoustic instability studied earlier by 



t from the nonlinear evolution of 



Malkov & Diamond 



the Drury's (|Malkov et al. 



(2009). Both 



120101) and Bells 's instabilities (after adding dissipation to the soliton 
solution) should result in shock-like nonlinear structures considerably shorter than the conven- 
tional CR precursor of the standard Bohm diffusion model. The magnetic field enhancement is 
clearly weaker in the acoustic case, as it is merely due to the individual shock compression in 
the instability generated shock-train. Besides, the soliton scenario is ex citing as it introduces 



these fascinating and ubiquitous objects (e.g., 



Ablowitz & Segur 



19811) to the SNR physics. 



However, the dominant instability should be selected on a case by case basis by treating the 
alternatives in a specific shock environment. 
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Fig. 1. — Dispersive properties of the two types of solitons: the short scale (lower branch) 
and the long scale (upper brunch). The soliton wave number K is shown in the units of kj = 
ItiJ/cBq. 
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Fig. 2. — 5 x -component of the soliton solution as a function of dimensionless coordinate s = 
in units of B max (B^ ax = SnpoC 2 ): B x /B max = ^fwcos(®) shown with the solid line 
and the amplitude envelope, ±y/w (dashed line). The soliton Mach number Ma = C/Ca = 3 
corresponds to the amplitude parameter a = 1 — M A 2 =8/9. 




Fig. 3. — Phase © as a function of coordinate s, shown with the solid line. The amplitude of the 
soliton is shown with the dashed line. 
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Fig. 4. — The ratio of m agnetic fiel d energy to ambient density po as a function of shock 



velocity V s a dopted fro m (IVink 



scaling from (Bell 



2008|) (points with error ba rs). The 



5lue dotted line is the 



2005), both also taken 



120041) . while the red dashed line is that of (|Volk et all 
from the Vink's compilation. An example of calculations, described in Sec|4]is shown with the 
green line. The field energy declines beyond V s ~ 2.7 ■ 10 4 to vanish at Vq ~ 9 ■ 10 4 km/s (not 
shown in the plot). At lower V s , where the shock acceleration is inefficient, this dependence 



breaks down and shou 



Drury & Voelk 



1981 



d transition to a low-efficiency acceleration regime in a bistable fashion 



Malkov 



19971) (schematically shown with the magenta dotted line). 



